%% to do:  X prior,  X gls;  financial,  lag level

clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SPEC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% pick survey
sur = 'FF'; % Two blue chip surveys - financial and economics. We may want to add the other survey. Bottom line - we don't have the economic survey. 

%%%%%%% financial variables -- take event out of dummy period
exev_dummy = 1; % Subtracts the 2-day window.  

%%%%%%% minimum obs for each forecaster (if inf -> balanced panel)
MinObs = 24; % forecasters come in and out of the survey. This matters because we do the sandwich matrix, and therefore we need SD for each forecaster. 
 % Double check the robustness to this. 
%%%%%%%% add lag
add_lag = 1; %%set 0 if not

%%%%%%%% add lag
add_laglevel = 1; %%set 0 if not

%%%%%%%% add constant/fixed effects
add_constant = 1; %%set 0 if not
add_FE = 1; % for each forecaster

%% std errors type
%% = 0: OLS
%% = 1: diagonal (heterosk across fcsters)
%% = 2: PC equal loading - this means common component in the error and then idiosyncratic components. factor model where each has a loading of 1. Then you estimate white. 
%% = 3: PC - this is the principal component. Allowing for different loadings. Might need more MinObs. 

sandwich = 2; 

%% do GLS
if sandwich
  Gls = 0;
else
  Gls = 0;
end

%% shrinkage
shrink = 0; %.5; We are shrinking the number of regressors. Stein estimator. 0 is OLS, the bigger it is, the more you shrink. 

%%%%%%%% regressions for individual forecasters
do_individual = 0;

%%%%%%% show figures
show_figures = 0;

%%%%%%% show figures
skip_regress = 0;

%%%%%%%% select event dummies
sel_separate = 0; % 


%% horizon
hor_list = 1:6;

%%%%%%%% choose time period for regression
start_r = []; %% if empty, start from beginning
start_r = [2008 06 01];
end_r = [];   %% if empty, end at the end of the sample


debg = 0;
disp([ 'debug : ' num2str(debg,0) ' -- press return to continue']);
%pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%set(0,'DefaultFigureVisible','off');  % all subsequent figures "off"
addpath('general');
addpath('BlueChip_051415');
gpath = [sur 'results'];
if ~exist(gpath,'dir'); mkdir(gpath); end;
filen = [gpath '/results.txt']; 
if ~debg 
  fileID = fopen(filen,'w');
else
  fileID = 1; %% if want to see results on screen
end


%% loading the contents of BC*_bp.mat (see readme file in this directory) 
eval(['load BC' sur '_bp_merged_dummy_051415;']); % Double Check this as we re-load the new data
eval(['load BL' sur '_dummy_051415;']);

%% create date strings and concatenate issue dates
eval(['dtnum = BC' sur '_bp.dtnum;']);
eval(['d_dates = datevec(BC' sur '_bp.dtnum);']);
d_dates = d_dates(:,1:3);

date_str = [ num2str(d_dates(:,1)) repmat('/',[size(d_dates,1) 1])  num2str(d_dates(:,2))  repmat('/',[size(d_dates,1) 1]) num2str(d_dates(:,3)) ];
dates = 10^4*d_dates(:,1)+10^2*d_dates(:,2)+d_dates(:,3);


%% TOTALLY TEMPORARY 

%Save = BCFF_bp.FederalFundsRate.cons.forecast;
%Save = [d_dates, Save];
%csvwrite('FFR_forecast.csv', Save);


%% create date strings and concatenate issue dates
temp = csvread('ColDates_051415.csv');
dtcol = temp(:,2);
coldates = temp(:,1);


%% size of panel / number of participants
N = eval(['length(BC' sur '_bp.partArray)']);
T = length(dates);

fprintf(fileID,'\n size of panel: number of participants: %1.0f, number of forecasts: %1.0f \n ',[N T]);


%% figure out when forecast horizon changes
eval(['cQ = BC' sur '_bp.cQ;']);
cQ1 = [ 0; cQ(1:end-1,1)];
chQ = (cQ(:,1) ~= cQ1); % 
chQN = NaN*chQ; 
chQN((chQ==0)) = 0;

%% same as chQ
eval(['shiftQ = BC' sur '_bp.shiftQ;']); % Dummy that is 1 if the forecast is shifted

%%%%%%%%%%%%%%%%%%%%%%%%% get QE and fwd guidance event dummies (D)
ev = csvread('event_codes.csv',1,1);

ev_dates = ev(:,1:3);
fwd_dummy = ev(:,4);
QE_announce_dummy = ev(:,5);
QE_dummy = ev(:,6);
twist_dummy = ev(:,7);
GDP_news_dummy = ev(:,8);
inf_news_dummy = ev(:,9);
QE1_dummy = ev(:,10);
QE2_dummy = ev(:,11);
QE3_dummy = ev(:,12);
GDP_change_good = ev(:,13);
GDP_change_bad = ev(:,14);

ev_names = 10^4*ev_dates(:,1)+10^2*ev_dates(:,2)+ev_dates(:,3);

ev_Dummy = zeros(T,size(ev_names,1));
ev_time = zeros(size(ev_names,1),1);
ev_chQ = zeros(size(ev_names,1),1);

if sel_separate == 1
  
  sel_ev = []; %% if 'inf' select all events, if [] none 
  %sel_ev = [20110809 20120125 20120313 20120425 20120913];  %[20110809];%[20120913];%[20120125];%[20120125];%[20121212];%[20120125];%[20121212];%; [20120125 20110809 20120913];%
%   sel_ev = [20110809  20120125   20120913];
  sel_ev = [20110809];
%   sel_ev = [20120913];
  
elseif ~sel_separate
  %% does not work w ~separate, will need to fix
  sel_ev = [(fwd_dummy==1) (QE_announce_dummy == 1) (QE_dummy == 1) (twist_dummy == 1) (GDP_news_dummy == 1) (GDP_news_dummy == 3) (inf_news_dummy == 1)]; % These are the dummies for forward guidance and quantitative easing.
  dum_names = {'fwdguidance' 'QE_announce' 'QE_continue' 'twist' 'GDP_bad' 'GDP_good' 'inf'};

end

 % The point of this is to see if the event takes place before or after a
 % given collection date
for kk = 1 : size(ev_dates,1) % This is looping over all of the FOMC statement dates 
  
  dtdiff = dtcol -datenum(ev_dates(kk,:));
  dtdiff((dtdiff < 0 )) = inf; 
  [mindx idx] = min(dtdiff);
  %index of closest value
  
  ev_Dummy(idx,kk) = 1; % This creates a dummy at the time of the FOMC statement
  ev_time(kk) = idx; % This indexes the time of the FOMC statement
  ev_chQ(kk) = chQ(idx); %This indexes is this is a point where there is a change in the forecast horizon
  %disp([ev_dates(kk,:)]);  disp([idx  coldates(idx-1) coldates(idx) dates(idx)]);   %pause
end

% Duall0 = zeros(T+1,5);
% eval(['Duall0(1:T,1) = BC' sur '_bp.fgQEnone;']);    Duall0(T+1,1) = 0;
% eval(['Duall0(1:T,2) = BC' sur '_bp.fgQE1;']);       Duall0(T+1,2) = 1;
% eval(['Duall0(1:T,3) = BC' sur '_bp.otherQE12;']);   Duall0(T+1,3) = 2;
% eval(['Duall0(1:T,4) = BC' sur '_bp.fgQE3;']);       Duall0(T+1,4) = 3;
% eval(['Duall0(1:T,5) = BC' sur '_bp.otherQEnone;']); Duall0(T+1,5) = 4;
% [ifd jfd] = find(Duall0(1:T,:));
% fd = sortrows([ifd jfd],1);

%%%%%%%%%%%%%%%%%%%%%%%%% construct BL macro surprises
nu_macro = 27; % number of macro variables (These include the consumer confidence indexes)
nu_fin = 47; % number of financial variables 
nu_other = 0; % This is non-financial ? Uncertainty index 

%% get names of BL series - this is meant to double check the names of the variables
% This is done manually 
%BLmacro_names = importdata('BLmacro_names.csv');
[NUM,TXT,RAW] = xlsread('BL_description.xls');
BL_desc = TXT(:,[1 2 end]);
BL_desc_macro = BL_desc(1:nu_macro,:);

%% construct bloomberg surprise = actual - forecast 
eval(['BLlist = fieldnames(BL' sur ');']);
BLmacro_list = cell(length(BL_desc_macro),1);
BLmacro_list = BLlist(7+1:7+nu_macro);

fn1 = 'forecastcol';
fn2 = 'actualcol';
fn3 = 'dtrelcol';

Act = NaN(T,length(BLmacro_list));
Fcst = NaN(T,length(BLmacro_list));
Diff = NaN(T,length(BLmacro_list));
Rel_Date = NaN(T,length(BLmacro_list));

%% fill zero for missing values (in btw)?
fill_zero = 1; % NOTE: the NaN data are months between the quarters

fprintf(fileID,'\n BLmacro_ data availability\n ');
ShrDummies_BLmacro_diff = NaN(length(BLmacro_list),1);

%% we are excluding GDP and filling it out w 2:4 (quarterly real GDP)
for j = 1:length(BLmacro_list)
  reg = BLmacro_list{j};
  eval([' fcst = getfield(BL' sur '.' reg ',fn1);']);
  eval([' act = getfield(BL' sur '.' reg ',fn2);']);
  eval([' release_date = getfield(BL' sur '.' reg ',fn3);']);

  reldates = datevec(release_date);
  Rel_Date(:,j) = 10^4*reldates(:,1)+10^2*reldates(:,2)+reldates(:,3);
  
  Act(:,j) = act;
  Fcst(:,j) = fcst;  
  Diff(:,j) = act-fcst;
  
  %% combine into one GDP news (first)
  if any( j == [2:4] )
    Act(~isnan(Act(:,j)),1) = Act(~isnan(Act(:,j)),j);
    Fcst(~isnan(Act(:,j)),1) = Fcst(~isnan(Act(:,j)),j);  
    Diff(~isnan(Act(:,j)),1) = Diff(~isnan(Act(:,j)),j);
  end
  
  ShrDummies_BLmacro_diff(j) = nanstd(Diff(:,j)); % Note that the NAN are referring to the MONTHS within the quarter
  %      fprintf(1,'\n  %1.0f  %2.20s  \n',[ j BLmacro_list{j}]);   format long g;  disp( [dates coldates Rel_Date(:,j) fcst act ]);  pause;
  
  iXj = find(~isnan(Diff(:,j))); % This just indexes the month that the quarterly data is listed for
  
  %fprintf(1,'  %1.0f  %2.20s  ,',[ j BLmacro_list{j}]);  
  fprintf(fileID,'\n %1.0f  %2.20s: beginning: ',[ j BLmacro_list{j}]);
  fprintf(fileID,'   %1.0f  ',dates(iXj(1)));
  fprintf(fileID,'  , end: ');
  fprintf(fileID,'   %1.0f  ',dates(iXj(end)));
  fprintf(fileID,'  , missing in btw (fraction):  %1.2f  ',1 - length(iXj)/(iXj(end)-iXj(1)+1));
  fprintf(fileID,'  , different missing in act and forecast :  %1.0f  ',(length(find(isnan(act)))-length(find(isnan(fcst))) ));
  
  
  if fill_zero & any(isnan(Diff(iXj(1):iXj(end),j))); 
    fprintf(fileID,' fill zero for missing values');
    iXjj = find(isnan(Diff(:,j)));
    iXjj = iXjj(find( ( iXjj > iXj(1) )&( iXjj < iXj(end) ) ));
    Diff(iXjj,j) = 0;

    %        disp([iXj(1) iXj(end) length(iXjj) any(isnan(Diff(iXj(1):iXj(end),j)))  ]);  pause
    
  end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%% construct BL financial surprises
%% get names of BL series
%BLfin_names = importdata('BLfin_names.csv');

%% construct bloomberg surprise = actual - forecast 
BL_desc_fin = BL_desc(nu_macro+1:nu_macro+nu_fin,:);
BLfin_list = cell(length(BL_desc_fin),1);
BLfin_list(1:end) = BLlist(7+nu_macro+1:7+nu_macro+nu_fin);

fn1 = 'actualcolWindow';
fn2 = 'actualcolEvents';
fn3 = 'actualcolWsubE';

Fin = NaN(T,length(BLfin_list),2);
EvFin = NaN(T,length(BLfin_list));

fprintf(fileID,'\n BLfin_ data availability\n ');
ShrDummies_BLfin = NaN(length(BLfin_list),1);

for j = 1:length(BLfin_list)
  reg = BLfin_list{j};
  eval([' fin = getfield(BL' sur '.' reg ',fn1);']);
  eval([' evfin = getfield(BL' sur '.' reg ',fn2);']);
  eval([' diffin = getfield(BL' sur '.' reg ',fn3);']);

  
  Fin(:,j,1) = fin;
  EvFin(:,j) = evfin;  
  Fin(:,j,2) = diffin; % This subtracts the 2-day window
  
  Fin(isinf(fin),j,1) = NaN;
  EvFin(isinf(evfin),j) = NaN;  
  Fin(isinf(diffin),j,2) = NaN;
  
    
  ShrDummies_BLfin(j,:) = nanstd(Fin(:,j,exev_dummy+1)); 
       % fprintf(1,'\n  %1.0f  %2.20s  \n',[ j BLfin_list{j}]);   format long g;  disp( [dates coldates fin evfin diffin]);  disp([nanstd(Fin(:,j,1)); nanstd(Fin(:,j,2))]); pause;
  
  iXj = find(~isnan(Fin(:,j,exev_dummy+1)));
  
  %fprintf(1,'  %1.0f  %2.20s  ,',[ j BLfin_list{j}]);  
  fprintf(fileID,'\n %1.0f  %2.20s: beginning: ',[ j BLfin_list{j}]);
  fprintf(fileID,'   %1.0f  ',dates(iXj(1)));
  fprintf(fileID,'  , end: ');
  fprintf(fileID,'   %1.0f  ',dates(iXj(end)));
  fprintf(fileID,'  , missing in btw (fraction):  %1.2f  ',1 - length(iXj)/(iXj(end)-iXj(1)+1));
   
end


%%%%%%%%%%%%%%%%%%%%%%%%% construct BL ``other'' surprises
BL_desc_oth = BL_desc(nu_macro+nu_fin+1:nu_macro+nu_fin+nu_other,1);
BLoth_list = cell(length(BL_desc_oth),1);
BLoth_list(1:nu_other) = BLlist(7+nu_macro+nu_fin+1:7+nu_macro+nu_fin+nu_other);

fn1 = 'actualcol';

Oth = NaN(T,length(BLoth_list));

fprintf(fileID,'\n BLoth_ data availability\n ');
ShrDummies_BLoth = NaN(length(BLoth_list),1);

for j = 1:length(BLoth_list)
  reg = BLoth_list{j};
  eval([' oth = getfield(BL' sur '.' reg ',fn1);']);
  Oth(:,j,1) = oth;
  Oth(isinf(oth),j,1) = NaN;
  ShrDummies_BLoth(j,:) = nanstd(Oth(:,j)); 
       %fprintf(1,'\n  %1.0f  %2.20s  \n',[ j BLoth_list{j}]);   format long g;  disp( [dates coldates oth]);  disp([nanstd(Oth(:,j)); nanstd(Oth(:,j))]); pause;
  
  iXj = find(~isnan(Oth(:,j)));
  
  %fprintf(1,'  %1.0f  %2.20s  ,',[ j BLoth_list{j}]);  
  fprintf(fileID,'\n %1.0f  %2.20s: beginning: ',[ j BLoth_list{j}]);
  fprintf(fileID,'   %1.0f  ',dates(iXj(1)));
  fprintf(fileID,'  , end: ');
  fprintf(fileID,'   %1.0f  ',dates(iXj(end)));
  fprintf(fileID,'  , missing in btw (fraction):  %1.2f  ',1 - length(iXj)/(iXj(end)-iXj(1)+1));
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% select regressors (X)
fprintf(fileID,'\n\n regressors: Event Dummies ' );
str_name = [ ];

%%%%%%%% select event dummies
if sel_separate == 1

  if isinf(sel_ev) % Note that these are the dates for the FOMC announcement of forward guidance
    X = ev_Dummy; % This does ALL dummies
    %         dum_names = {num2str(ev_names)};
    dum_names = cell(length(ev_names),1);
    dum_ev_chQ = zeros(length(ev_names),1);
    for kk = 1 : length(ev_names)
      dum_names(kk) = {num2str(ev_names(kk))};
      dum_ev_chQ(kk) = ev_chQ(kk);
    end
     
  elseif ~isempty(sel_ev)
    %% put them in ascending order
    
    X = zeros(T,length(sel_ev));
    dum_ev_chQ = zeros(length(sel_ev),1);
    for kk = 1 : length(sel_ev)
      X(:,kk) = ev_Dummy(:,find(ev_names == sel_ev(kk)));
      dum_names(kk) = {num2str(sel_ev(kk))};
      dum_ev_chQ(kk) = ev_chQ(find(ev_names == sel_ev(kk)));
    end
  
  elseif isempty(sel_ev)
    X = [];
    dum_names = {};
    dum_ev_chQ = [];
  end  
  
elseif ~sel_separate
  if size(sel_ev,2) ~= length(dum_names)
    error('size(sel_ev,2) ~= length(dum_names)');
  end  
  X = zeros(T,size(sel_ev,2));
  for kk = 1 : size(sel_ev,2)        
    X(:,kk) = sum(ev_Dummy(:,find(sel_ev(:,kk))),2); 
%    dum_ev_ChQ(kk) =  
  end
  
end
%% make sure there are no duplicate dummies
if  ~isempty(sel_ev)
  if rank(X) < size(X,2)  
    [iid jjd] = find(X);
    jjd = jjd(find(diff(iid) == 0)+1);
    for kk2= 1:length(jjd)
      shft = 1;
      while any(jjd == jjd(kk2)-shft)
        shft = shft+1;
      end      
      dum_names(jjd(kk2)-shft) = { [dum_names{jjd(kk2)-shft} '-' dum_names{jjd(kk2)}] };
    end
    X(:,jjd) = [];
    dum_names(jjd) = [];
    dum_ev_chQ(jjd) = [];
  end 
end


if ~isinf(sel_ev)
  for kk = 1 : size(X,2)
    fprintf(fileID,' %2.19s , ', dum_names{kk});
    str_name = [str_name '_' dum_names{kk}];
  end

elseif isinf(sel_ev)
  fprintf(fileID,' %2.19s , All Dates ');
  str_name = [str_name '_AllDates' ];
end


ev_loc = size(X,2);
ShrDummies = zeros(ev_loc,1);



%%%%%%%% select macro news
%%  1  RGDP_First  ,  2  RGDP_Second  ,  3  RGDP_Third  ,  4  PCEDEFM  ,  5  PCECMOM  ,  6  GDPCTOT  ,  7  PITL  ,  8  NAPMPMI  ,  9  NAPMNMI  ,  10  CNSTTMOM  ,  11  USTBTOT  ,  12  NFP  ,  13  USURTOT  ,  14  RSTAMOM  ,  15  RSTAXMOM  ,  16  CPI  ,  17  CPTI  ,  18  IP  ,  19  CPUP  ,  20  DGNO , 21 CONCCONF , 22 INJCJC , 23 LEI , 24  NHSPSTOT , 25  NHSLTOT, 26  PPI , 27 CONSSENT
BLmacro__sel = [ ];
%BLmacro__sel = [2:4 8:9 12:21];
%BLmacro__sel = [3 8 12:13 15:20]; BASELNE
BLmacro__sel = [17 21 16 1 22 23 8 25 12 14 13]; %GSS Variables
fprintf(fileID,'BLmacro_ regressors: ' );

BLmacro__names = cell(length(BLmacro__sel),1);
for kk = 1 : length(BLmacro__sel)
  BLmacro__names(kk) = {BLmacro_list{BLmacro__sel(kk)}};
  fprintf(fileID,' %2.19s  ', BLmacro__names{kk});
    fprintf(fileID,' & %2.999s  ', BL_desc_macro{BLmacro__sel(kk),2});
        fprintf(fileID,' & %2.999s  \\\\ \n ', BL_desc_macro{BLmacro__sel(kk),3});
  str_name = [str_name '_' BLmacro__names{kk}];
end


%%%%%%%% select fin news
%% 1  BBOX  ,  2  GNMA15  ,  3  GNMA30  ,  4  FMNMA15  ,  5  FMNMA30  ,  6  GOLD15  ,  7  GOLD30  ,  8  MBS15  ,  9  MBS30  ,  10  FNMGVN5  ,  11  FNMGVN10  ,  12  FNMGVN30  ,  13  SP500  ,  14  DJIA  ,  15  FX_EDUS  ,  16  BE1  ,  17  BE5  ,  18  BE10  ,  19  BE20  ,  20 BE30  ,  21  DGS1  ,  22  DGS5  ,  23  DGS10  ,  24  DGS20  ,  25  DGS30  ,  26  TIPS10  ,  27  TIPS30  ,  28  MoodysAAA10  ,  29  MoodysBAA10  ,  30  VIX 
%BLfin__sel = [1 9 13 15 23 26];
BLfin__sel = [9 13 15 23 26];
%BLfin__sel = [ ];

fprintf(fileID,'BLfin_ regressors: ' );

BLfin__names = cell(length(BLfin__sel),1);
for kk = 1 : length(BLfin__sel)
  BLfin__names(kk) = {BLfin_list{BLfin__sel(kk)}};
  fprintf(fileID,' %2.19s  ', BLfin__names{kk});
    fprintf(fileID,' & %2.999s  ', BL_desc_fin{BLfin__sel(kk),2});
    fprintf(fileID,' & %2.999s  \\\\ \n ', BL_desc_fin{BLfin__sel(kk),3});

  str_name = [str_name '_' BLfin__names{kk}];
end


%%%%%%%% select oth news
%1  SEPUI 
BLoth__sel = [];
%BLoth__sel = [ ];

fprintf(fileID,'BLoth_ regressors: ' );

BLoth__names = cell(length(BLoth__sel),1);
for kk = 1 : length(BLoth__sel)
  BLoth__names(kk) = {BLoth_list{BLoth__sel(kk)}};
  fprintf(fileID,' %2.19s  ', BLoth__names{kk});
    fprintf(fileID,' & %2.999s  ', BL_desc_oth{BLoth__sel(kk),2});
    fprintf(fileID,' & %2.999s \\\\ \n', BL_desc_oth{BLoth__sel(kk),3});

    str_name = [str_name '_' BLoth__names{kk}];
end

%% combine the covariates - note that before we put it all together, X is the event dummy
X = [X Diff(:,BLmacro__sel) Fin(:,BLfin__sel,exev_dummy+1)  Fin(:,BLoth__sel)];
ShrDummies = [ShrDummies; ShrDummies_BLmacro_diff(BLmacro__sel); ShrDummies_BLoth(BLoth__sel)]; 

if shrink
  f_sh = find(ShrDummies);
  X_d0 = zeros(length(f_sh),size(X,2));
  X_d0(:,f_sh) = diag(shrink*ShrDummies(f_sh));
else
  X_d0 = 0;
end
        
reg_names = {dum_names{:} BLmacro__names{:} BLfin__names{:}  BLoth__names{:}};


%%%%%%%% modify name with variosu options
if add_lag; 
  fprintf(fileID,'  lag ');
  str_name = [str_name '_lag']; 
  reg_names = {reg_names{:} 'lag'};
end;

if add_laglevel; 
  fprintf(fileID,'  lag level');
  str_name = [str_name '_laglev']; 
  reg_names = {reg_names{:} 'laglev'};
end;

if add_constant; 
  fprintf(fileID,'  constant ');
  str_name = [str_name '_const']; 
  if add_FE
    fprintf(fileID,'  Fixed Effects ');
    str_name = [str_name '_FE']; 
  end
end;

if sandwich
  fprintf(fileID,'  Sandwich St Errors');
  str_name = [str_name '_SWICH']; 

  if sandwich == 1
    fprintf(fileID,'  diagonal ');
    str_name = [str_name '_diag']; 
    
  elseif sandwich == 2
    fprintf(fileID,'  principal component ');
    str_name = [str_name '_PC'];     
  end
  
  if Gls
    fprintf(fileID,'  GLS ');
    str_name = [str_name '_GLS'];     
  end
    
end

if shrink
  fprintf(fileID,'  shrinkage %2.2f',shrink);
  str_name = [str_name '_shrink' num2str(10*shrink,'%2.0f')]; 
end


%%%%%%%% choose time period for regression
if ~isempty(start_r)
  i_s = find(sum(abs(d_dates(:,1:2) - repmat(start_r(1:2),[T 1])),2) == 0);
  str_name = [str_name '_START' num2str(dates(i_s))];
  fprintf(fileID,'\n\n Beginning of sample:   %2.0f , ',  dates(i_s));
else
  i_s = 1;
end

if ~isempty(end_r)
    i_e = find(sum(abs(d_dates(:,1:2) - repmat(end_r(1:2),[T 1])),2) == 0);
    str_name = [str_name '_START' num2str(dates(i_e))];
    fprintf(fileID,'  End of sample:   %2.0f \n\n ',  dates(i_e));
else
    i_e = T;
end

%str_name = [str_name '_MinObs' num2str(MinObs)];
fprintf(fileID,'  Minimum number of observations for each forecaster:   %2.0f \n\n ',  MinObs);
fprintf(fileID,'  \n\n %2.999s \n\n ',  str_name);


X = X(i_s:i_e,:);    

%% adjust event dates
ev_iijj = find((ev_time >= i_s) & (ev_time <= i_e));
ev_dates = ev_dates(ev_iijj,:);
ev_names = ev_names(ev_iijj);
ev_time = ev_time(ev_iijj)-i_s+1;

% Adjusting dates for the plot 

%% create directory with regression name
if ~debg
  gpath_r = [gpath str_name];
  if ~exist(gpath_r,'dir'); mkdir(gpath_r); end;
  filenn = [gpath_r '/results' str_name '.txt'];
  eval(['!cp ' filen '  ' filenn]);
  fileID = fopen(filenn,'a');
elseif debg
  gpath_r = gpath;
end  

if  do_individual
    
  gpath_r_ind = [gpath_r '/ind'];
  if ~exist(gpath_r_ind,'dir'); mkdir(gpath_r_ind); end;
  filen_ind = [ gpath_r_ind '/results_ind.txt'];   
  eval(['!cp ' filen '  ' filen_ind]);
  fileID_ind = fopen(filen_ind,'a');

end

gpath_r0 = gpath_r; 
if ~exist(gpath_r,'dir'); mkdir(gpath_r); end;

%% identify/eliminate periods for which we do not have obs for the regressor
iX = ( sum(isnan(X),2) > 0 );
tT = size(X,1)-sum(iX);
X(iX,:) = [];

fprintf(fileID,'\n\n regression period, number time periods: %1.0f, beginning: ',[tT]);
fprintf(fileID,'   %1.0f  ',dates(i_s));
fprintf(fileID,'  , end: ');
fprintf(fileID,'   %1.0f  ',dates(i_e));


%%%%%%%%%%%%%%%%%%% pick forecast (Y)
eval(['fcstlist = fieldnames(BC' sur '_bp);']);
fcstlist = fcstlist(3:end-9);
flist = [17:19 1 5 9 10]; 

bounds = [ -2,2; NaN,NaN;NaN,NaN;NaN,NaN;-1,0.5;NaN,NaN;NaN,NaN;NaN,NaN;-1.6,1;-1.6,1;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;-2,2;-2,2;-1,2.5];  
bounds_policy = [ -2,2; NaN,NaN;NaN,NaN;NaN,NaN;-0.5,0.5;NaN,NaN;NaN,NaN;NaN,NaN;-0.5,0.5;-0.5,0.5;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;NaN,NaN;-1,1;-2,2;-1,1.2];  


for  v = flist; % This loops over forecast variables 
  
  fcst = fcstlist{v}; 
  %fcst = 'GDP'

  for kk = 1:ev_loc
    eval(['ev_effect_' dum_names{kk}  '.' fcst  '= NaN(length(hor_list)+1,2,2);']);
  end
  fprintf(fileID,'\n\n\n VARIABLE BEING FORECASTED: %2.20s \n ',fcst);

  fn = 'forecast'; 
  %% consensus
  eval([' y_cons = getfield(BC' sur '_bp.' fcst '.cons,fn);']);

  Y_panel = zeros(T,N,6); % Constructing a panel of forecasters
  for i = 1:N 
    %% each forecaster
    eval(['Y_panel(:,i,:) = getfield(BC' sur '_bp.' fcst '.f_' num2str(i) ',fn);']);
  end
  
  for hor = hor_list 

    % disp([y_cons(:,hor) nanmean(Y_panel(:,:,hor),2)]);     pause

    %% y lagged: we use same quarter when chQW(t)=0, and next quarter when  chQW(t)=1 (NaN for hor =6 -- we do not have the 'next quarter')
    if hor < 6 
      y_lag = (1-chQ).*[NaN; y_cons(1:end-1,hor)] + chQ.*[NaN; y_cons(1:end-1,hor+1)];
      Y_lag = repmat((1-chQ),[1 N]).*[NaN*ones(1,N); Y_panel(1:end-1,:,hor)] +  repmat(chQ,[1 N]).*[NaN*ones(1,N); Y_panel(1:end-1,:,hor+1)];
    elseif hor == 6
      y_lag = (1-chQ).*[NaN; y_cons(1:end-1,hor)] + chQ.*chQN;
      Y_lag = repmat((1-chQ),[1 N]).*[NaN*ones(1,N); Y_panel(1:end-1,:,hor)] +  repmat(chQ,[1 N]).*repmat(chQN,[1 N]);
    end
    
    y_diff = y_cons(:,hor) - y_lag; % Note that here we are using the change in the forecast (either the concesus or the individual)
    Y_diff = Y_panel(:,:,hor) - Y_lag;
    
    %% chop off to chosen sample
    y_diff = y_diff(i_s:i_e,:);
    Y_diff = Y_diff(i_s:i_e,:);
    y_lag = y_lag(i_s:i_e,:);
    Y_lag = Y_lag(i_s:i_e,:);
                                      
    fprintf(fileID,'\n\n HORIZON: %2.0f,  consensus, difference \n',hor);
    
    %% consensus
    Y_reg = y_diff;
    X_reg = X;
    Resid_ev = NaN*ones(size(Y_reg));

    % eliminate periods for which X is NaN
    Y_reg(iX,:) = [];
    
    X_d = X_d0;
    if add_lag
      y_diff1 = [NaN; y_diff(1:end-1)];  
      y_diff1(iX,:) = [];
      X_reg = [X_reg , y_diff1];

      if shrink
        X_d = [X_d zeros(size(X_d,1),size(y_diff1,2))];
      end
    end  
    if add_laglevel 
      y_lag1 = y_lag;  
      y_lag1(iX,:) = [];
      X_reg = [X_reg , y_lag1];

      if shrink
        X_d = [X_d zeros(size(X_d,1),size(y_lag1,2))];
      end
    end      
    if add_constant
      X_reg = [X_reg , ones(size(X_reg,1),1)];
      if shrink
        X_d = [X_d zeros(size(X_d,1),1)];
      end
    end    
    reg_names_reg = reg_names;
    
    fileID0 = fileID;
    ispanel = 0;
    run_ols;
    if sel_separate 
        for kk = 1:ev_loc
          if strcmp(dum_names{kk},reg_names_reg{kk})
            eval(['ev_effect_' dum_names{kk}  '.' fcst  '(hor+dum_ev_chQ(kk),:,1) = [beta(kk),stdev(kk)];']);
          end
        end
    end
    if ~sel_separate 
        for kk = 1:ev_loc
          if strcmp(dum_names{kk},reg_names_reg{kk})
            eval(['ev_effect_' dum_names{kk}  '.' fcst  '(hor,:,1) = [beta(kk),stdev(kk)];']);
          end
        end
    end
    Resid_ev(~iX,:) = Resid_ev0;
    

    
    %% panel
    Y_reg = Y_diff;
    Resid_ev = NaN*ones(size(Y_reg));
    
    % eliminate periods for which X is NaN
    Y_reg(iX,:) = [];
    X_reg = X;
    % eliminate forecasters for whom we  less than 3 observations 
    %jY = all(isnan(Y_reg));
    minobs = min(MinObs,size(Y_reg,1));
    jY = (size(Y_reg,1)-sum(isnan(Y_reg))<=minobs);
    Y_reg(:,jY) = [];

    nN = size(Y_reg,2);
    X_d = nN*shrink*X_d0;
    
    fprintf(fileID,'\n\n HORIZON: %2.0f,  panel (%2.0f forecasters, %2.0f periods), difference \n',[hor nN tT]);
    X_reg = kron(ones(nN,1),X_reg);
    if add_lag

      Y_diff1 = [NaN(1,size(Y_diff,2)); Y_diff(1:end-1,:)];  
      Y_diff1(iX,:) = [];
      Y_diff1(:,jY) = [];
      
      add_X = Y_diff1(:);
      % if want to have separate coefs
      %add_X = kron(eye(nN),ones(tT,1));
      %add_X(find(add_X)) = Y_diff1;
      X_reg = [X_reg, add_X];     
      %% disp('check'); for jjj=1:nN; disp([add_X(1+(jjj-1)*tT:jjj*tT,jjj) Y_diff1(:,jjj)]); disp(any(add_X(1+(jjj-1)*tT:jjj*tT,jjj)-Y_diff1(:,jjj))); pause; end;
      if shrink
        X_d = [X_d zeros(size(X_d,1),size(add_X,2))];
      end
    end
    if add_laglevel

      Y_lag1 = Y_lag;  
      Y_lag1(iX,:) = [];
      Y_lag1(:,jY) = [];
      
      add_X = Y_lag1(:);
      X_reg = [X_reg, add_X];     
      %% disp('check'); for jjj=1:nN; disp([add_X(1+(jjj-1)*tT:jjj*tT,jjj) Y_lag1(:,jjj)]); disp(any(add_X(1+(jjj-1)*tT:jjj*tT,jjj)-Y_lag1(:,jjj))); pause; end;
      if shrink
        X_d = [X_d zeros(size(X_d,1),size(add_X,2))];
      end
    end
    if add_constant 
      if add_FE
        add_X = kron(eye(nN),ones(tT,1));
      else
        add_X = ones(tT*nN,1);      
      end
      X_reg = [X_reg, add_X];      
      if shrink
        X_d = [X_d zeros(size(X_d,1),size(add_X,2))];
      end      
    end
    Y_reg = Y_reg(:);
    reg_names_reg = reg_names;

    fileID0 = fileID;
    ispanel = 1;
    run_ols;
    
    if sel_separate 
        for kk = 1:ev_loc
          if strcmp(dum_names{kk},reg_names_reg{kk})
            eval(['ev_effect_' dum_names{kk}  '.' fcst  '(hor+dum_ev_chQ(kk),:,2) = [beta(kk),stdev(kk)];']);
          end
        end
    end
    if ~sel_separate 
        for kk = 1:ev_loc
          if strcmp(dum_names{kk},reg_names_reg{kk})
            eval(['ev_effect_' dum_names{kk}  '.' fcst  '(hor,:,2) = [beta(kk),stdev(kk)];']);
          end
        end
    end    
    
    Resid_ev(~iX,~jY) = reshape(Resid_ev0,[tT nN]);
    if nN > 0 
      Std0 = NaN(nN,1);
      Std0(~jY) = nanvar(Resid0,1).^.5;
    else
      Std0 = NaN;
    end
    
    %% Individual forecasters
    if do_individual

      fprintf(fileID_ind,'\n\n\n VARIABLE BEING FORECASTED: %2.20s \n ',fcst);
      fprintf(fileID_ind,'\n\n HORIZON: %2.0f \n  individual forecasters  \n',hor);
      beta_all = NaN(size(X,2)+add_lag,size(Y_diff,2));
      stdev_all = NaN(size(X,2)+add_lag,size(Y_diff,2));
      info_all = NaN(size(Y_diff,2),4);
      
      for jj = 1:size(Y_diff,2)
        Y_reg = Y_diff(:,jj);
        % eliminate periods for which X is NaN
        Y_reg(iX,:) = [];
        % eliminate forecasters for whom we have fewer than 'excl' observations 
        if ~(size(Y_reg,1)-sum(isnan(Y_reg))< excl) %%~all(isnan(Y_reg));

          fprintf(fileID_ind,'\n %2.0f: ',jj);  
          X_reg = X;
          X_d = X_d0;
          if add_lag
            y_diff1 = [NaN; Y_diff(1:end-1,jj)];  
            y_diff1(iX,:) = [];
            X_reg = [X_reg , y_diff1];        
            if shrink
              X_d = [X_d zeros(size(X_d,1),size(y_diff1,2))];
            end
          end  
          if add_laglevel 
            y_lag1 = Y_lag(:,jj);  
            y_lag1(iX,:) = [];
            X_reg = [X_reg , y_lag1];
            if shrink
              X_d = [X_d zeros(size(X_d,1),size(y_lag1,2))];
            end
          end     
          if add_constant
            X_reg = [X_reg , ones(size(X_reg,1),1)];
            if shrink
              X_d = [X_d zeros(size(X_d,1),1)];
            end
          end 
          reg_names_reg = reg_names;
          fileID0 = fileID_ind;
          ispanel = 0;
          run_ols; 
          beta_all(:,jj) = beta(1:size(X,2)+add_lag);
          stdev_all(:,jj) = stdev(1:size(X,2)+add_lag);
          info_all(jj,1) = nanstd(Y_diff(:,jj));
          %[jj info_all(jj,1)]
          info_all(jj,2:end) = info;
        end
        
      end
      beta_all = nanmean(beta_all,2);
      stdev_all = nanmean(stdev_all,2);
      fprintf(fileID,'\n\n  individual forecasters  \n',hor);
      for jj = 1:length(reg_names)
        fprintf(fileID,' %2.9s: ',[reg_names{jj}]);
        fprintf(fileID,' %2.4g (%2.4g) ; ',[beta_all(jj) stdev_all(jj)  ]);
      end
      
    end
    
    if show_figures
      for kk = 1 : size(ev_dates,1)      
        
        yy = Y_diff(ev_time(kk),:);
        rr = Resid_ev(ev_time(kk),:);
        fig_n = 100+kk;
        plot_bar;     
        
        fig_n = 200+kk;
        if ~all(isnan(rr))
          figure(fig_n)
          subplot(3,2,hor);
          P = plot(yy,rr,'x');
          %hold on;        
          %P1 = plot([nanmean(yy),nanmean(yy)],get(gca,'YLim'),'-.');
          hold on;
          P2 = plot(get(gca,'XLim'),[nanmean(rr),nanmean(rr)],'-.');
          hold off;
          title(['horizon: ' num2str(hor,' %2.0f')] );
        end
        
        yy = Std0;
        rr = Resid_ev(ev_time(kk),:);
        fig_n = 300+kk;
        if ~all(isnan(rr))
          figure(fig_n)
          subplot(3,2,hor);
          P = plot(yy,rr,'x');
          %hold on;        
          %P1 = plot([nanmean(yy),nanmean(yy)],get(gca,'YLim'),'-.');
          hold on;
          P2 = plot(get(gca,'XLim'),[nanmean(rr),nanmean(rr)],'-.');
          hold off;
          title(['horizon: ' num2str(hor,' %2.0f')] );
        end

        yy = Resid_ev(ev_time(kk),:)./Std0';
        rr = Resid_ev(ev_time(kk),:);
        fig_n = 400+kk;
        if ~all(isnan(rr))
          figure(fig_n)
          subplot(3,2,hor);
          P = plot(yy,rr,'x');
          %hold on;        
          %P1 = plot([nanmean(yy),nanmean(yy)],get(gca,'YLim'),'-.');
          hold on;
          P2 = plot(get(gca,'XLim'),[nanmean(rr),nanmean(rr)],'-.');
          hold off;
          title(['horizon: ' num2str(hor,' %2.0f')] );
        end

        
        if do_individual
          yyy = info_all(:,1);
          rrr = info_all(:,2);
          fig_n = 1001;   
          figure(fig_n)
          subplot(3,2,hor); 
          P = plot( info_all(:,1), info_all(:,2),'x'); 
          title(['horizon: ' num2str(hor,' %2.0f')] );
          
          fig_n = 1002;      
          figure(fig_n)
          subplot(3,2,hor); 
          P = plot( info_all(:,4), info_all(:,3),'x'); 
          title(['horizon: ' num2str(hor,' %2.0f')] );
          
          fig_n = 1003;      
          figure(fig_n)
          subplot(3,2,hor); 
          P = plot( info_all(:,1), info_all(:,3),'x'); 
          title(['horizon: ' num2str(hor,' %2.0f')] );

        end
        
      end 
 
      
    end %% end of showfigure
    
%       resid_average = nanvar(Resid_ev,1,2); % variance across forecasters
%       dates_label = dates(i_s:i_e);
%       dates_label(iX) =[];
%       resid_average(iX) = [];
%       ind = isnan(resid_average);
%       resid_average(ind) = [];
%       dates_label(ind) = [];
%       dates_label = round(dates_label/10000);
%       fig_n = 1;      
%       figure(fig_n)
%       P = plot(1:length(resid_average), resid_average);
%       set(P,'LineWidth',2,'Color',[0 0 0],'LineStyle','-');
%       hold on;
%       
%       set(gca,'XTick',1:8:length(resid_average));
%       set(gca,'XTickLabel',num2str(dates_label(1:8:end)), 'fontsize',18);
%       saveas(figure(fig_n),[gpath_r0 '/residuals_hor' num2str(hor) '_var_' fcstlist{v} ],'eps')
%       close(figure(fig_n))     
    
  end %% end horizon 

  for kk = 1:ev_loc
%    load([ gpath_r0 '/ev_effect_' dum_names{kk} '_' fcst '.mat '],'-mat')
    save([ gpath_r0 '/ev_effect_' dum_names{kk} '_' fcst '.mat '],['ev_effect_' dum_names{kk}]) 
  end
  
  for kk = 1:ev_loc
    nm_g = {'con' 'pan'};
    for gg = 1:2 %% 1 is cons, 2 is panel
      % mean irf      
      eval(['m_irf = ev_effect_' dum_names{kk}  '.' fcst ' (:,1,gg); ']);
      % stdev irf        
      eval(['std_irf = ev_effect_' dum_names{kk}  '.' fcst ' (:,2,gg); ']);
      fig_n = kk;      
      figure(fig_n)
      P = plot((0:length(hor_list)), m_irf);
      set(P,'LineWidth',4,'Color',[0 0 0],'LineStyle','-');
      hold on;
      P = plot((0:length(hor_list)), m_irf+.994*std_irf,(0:length(hor_list)), m_irf-.994*std_irf);
      set(P(:),'LineWidth',3,'Color',[0 0 0],'LineStyle','-.');
      hold on;
      P = plot((0:length(hor_list)), m_irf+1.96*std_irf,(0:length(hor_list)), m_irf-1.96*std_irf);
      set(P(:),'LineWidth',3,'Color',[0 0 0],'LineStyle',':');
      hold off;
      if sel_separate
    	 eval(['axis([0 6 ' num2str(bounds(v,1)) ' ' num2str(bounds(v,2)) ']);']);
      end
      if ~sel_separate
         eval(['axis([0 4 ' num2str(bounds_policy(v,1)) ' ' num2str(bounds_policy(v,2)) ']);']);
      end
      set(gca,'XTick',0:1:length(hor_list)+1);
      set(gca,'XTickLabel',num2str((0:length(hor_list))'), 'fontsize',18);

      saveas(figure(fig_n),[gpath_r0 '/ev_effect_' dum_names{kk} '_' fcst '_' nm_g{gg}],'eps')
      close(figure(fig_n))
      
    end
    
  end
  
  % Caluclating and plotting the average residual across forecasters over time
        

%     if ~sel_separate
%         for gg = 1:2 %% 1 is cons, 2 is panel        
%           for kk = 1:ev_loc 
%               eval(['m_irf_' dum_names{kk} '= ev_effect_' dum_names{kk}  '.' fcst ' (:,1,gg); ']);      
%               eval(['std_irf_' dum_names{kk} '= ev_effect_' dum_names{kk}  '.' fcst ' (:,2,gg); ']);
%           end
% 
%           % Manually adding up the effects for now - NEED TO AUTOMATE THIS
% 
% %           std_irf_20110809 = std_irf_fwdguidance;
% %           m_irf_20120125 = m_irf_fwdguidance + m_irf_QE2;
% %           std_irf_20120125 = std_irf_fwdguidance + std_irf_QE2;
% %           m_irf_20120913 = m_irf_fwdguidance + m_irf_QE3;
% %           std_irf_20120913 = std_irf_fwdguidance + std_irf_QE3;
%           
%           for xx = 1:length(fwg_events)
%               eval(['m_irf_' fwg_events{xx} '= m_irf_fwdguidance*' num2str(fwg_events_ind(1,xx)) '+ m_irf_QE_announce*' num2str(fwg_events_ind(2,xx))  '+ m_irf_QE_continue*' num2str(fwg_events_ind(3,xx))  '+ m_irf_twist*' num2str(fwg_events_ind(4,xx))  '+ m_irf_GDP_bad*' num2str(fwg_events_ind(5,xx))  '+ m_irf_GDP_medium*' num2str(fwg_events_ind(6,xx)) '+ m_irf_GDP_good*' num2str(fwg_events_ind(7,xx)) '+ m_irf_inf*' num2str(fwg_events_ind(8,xx)) ';']);              
%               %eval(['std_irf_' fwg_events{xx} '= std_irf_fwdguidance*' num2str(fwg_events_ind(1,xx)) '+ std_irf_QE_announce*' num2str(fwg_events_ind(2,xx))  '+ std_irf_QE_continue*' num2str(fwg_events_ind(3,xx))  '+ std_irf_twist_dummy*' num2str(fwg_events_ind(4,xx))  '+ std_irf_GDP_news_dummy*' num2str(fwg_events_ind(5,xx))  '+ std_irf_inf_news_dummy*' num2str(fwg_events_ind(6,xx)) ';']);              
%               fig_n = xx;      
%               figure(fig_n)
%               eval(['P = plot((0:length(hor_list)), m_irf_' fwg_events{xx} ')' ]);
%               set(P,'LineWidth',2,'Color',[0 0 0],'LineStyle','-');
%               hold on;
%               %eval(['P = plot((0:length(hor_list)), m_irf_' fwg_events{xx} '+.994*std_irf_' fwg_events{xx} ',(0:length(hor_list)), m_irf_' fwg_events{xx} '-.994*std_irf_' fwg_events{xx} ')']);
%               %set(P(:),'LineWidth',1,'Color',[0 0 0],'LineStyle','-.');
%               %hold on;
%               %eval(['P = plot((0:length(hor_list)), m_irf_' fwg_events{xx} '+1.96*std_irf_' fwg_events{xx} ',(0:length(hor_list)), m_irf_' fwg_events{xx} '-1.96*std_irf_' fwg_events{xx} ')']);
%               %set(P(:),'LineWidth',1,'Color',[0 0 0],'LineStyle',':');
%               %hold off;
%               eval(['axis([0 6 ' num2str(bounds(v,1)) ' ' num2str(bounds(v,2)) ']);']); 
%               set(gca,'XTick',0:1:length(hor_list)+1);
%               set(gca,'XTickLabel',num2str((0:length(hor_list))'));  
%               
%               saveas(figure(fig_n),[gpath_r0 '/ev_effect_added' fwg_events{xx} '_' fcst '_' nm_g{gg}],'eps')
%               close(figure(fig_n))
%           end
% 
%         end 
%     end

  if show_figures
    for kk = 1 : size(ev_dates,1)
      
      figure(100+kk)
      subtitle([ num2str(ev_names(kk),' %2.0f') ' - ' fcst]);
      saveas(figure(100+kk),[gpath_r '/FigHist' num2str(ev_names(kk),' %2.0f') '_' fcst],'eps')
      close(figure(100+kk))
      
      figure(200+kk)
      subtitle([ num2str(ev_names(kk),' %2.0f') ' - ' fcst]);
      saveas(figure(200+kk),[gpath_r '/FigScatt' num2str(ev_names(kk),' %2.0f') '_' fcst],'eps')
      close(figure(200+kk))
      
      figure(300+kk)
      subtitle([ num2str(ev_names(kk),' %2.0f') ' - ' fcst]);
      saveas(figure(300+kk),[gpath_r '/Fig-VScatt' num2str(ev_names(kk),' %2.0f') '_' fcst],'eps')
      close(figure(300+kk))
      
      
      figure(400+kk)
      subtitle([ num2str(ev_names(kk),' %2.0f') ' - ' fcst]);
      saveas(figure(400+kk),[gpath_r '/Fig-TstScatt' num2str(ev_names(kk),' %2.0f') '_' fcst],'eps')
      close(figure(400+kk))
      
    end

    if do_individual
      figure(1001)
      subt = [ 'var_vs_s2-' fcst];
      subtitle(subt);    
      saveas(figure(1001),[gpath_r_ind '/Fig' subt ],'eps')
      
      figure(1002)
      subt = [ 'nobs_vs_s2overvar-' fcst];
      subtitle(subt);    
      saveas(figure(1002),[gpath_r_ind '/Fig' subt ],'eps')
      
      figure(1003)
      subt = [ 'var_vs_s2overvar-' fcst];
      subtitle(subt);    
      saveas(figure(1003),[gpath_r_ind '/Fig' subt ],'eps')
    end    
    close all; 
  end
end  %% end fcst variable



